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Abstract. We present a simple non-equilibrium model of mass condensation with 
Lennard-Jones interactions between particles and the substrate. We show that when 
some number of particles is deposited onto the surface and the system is left to 
equilibrate, particles condense into an island if the density of particles becomes higher 
than some critical density. We illustrate this with numerically obtained phase diagrams 
for three-dimensional systems. We also solve a two-dimensional counterpart of this 
model analytically and show that not only the phase diagram but also the shape 
of the cross-sections of three-dimensional condensates qualitatively matches the two- 
dimensional predictions. Lastly, we show that when particles are being deposited with 
a constant rate, the system has two phases: a single condensate for low deposition 
rates, and multiple condensates for fast deposition. The behaviour of our model is thus 
similar to that of thin film growth processes, and in particular to Stranski-Krastanov 
growth. 
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1. Introduction 

Non-equilibrium statistical mechanics has witnessed a rapid progress in recent years, and 
has been applied to a variety of problems in physics, chemistry, biology, economy, and 
social sciences. However, in contrast to equilibrium systems, which can be conveniently 
studied by using the concept of the statistical ensemble, a unihed theoretical framework 
applicable to all non-equilibrium systems does not exist, and whether such a framework 
will eventually emerge remains to be seen. 

Despite that, signihcant progress has been made in the last two decades for a class of 
models called driven diffusive systems [I] which - even though being far from equilibrium 
- can be studied within the same statistical ensemble framework as equilibrium models. 
These models share a common feature: the steady-state probability of a microstate can 
be expressed analytically as a function of transition rates which define the dynamics 
of the model. Examples of such systems are the zero-range process (ZRP) [21 El 0], 
closely related to its equilibrium counterpart: balls-in-boxes model (B-in-B) [S], the 
asymmetric simple exclusion process (ASEP) |6] and its totally asymmetric version 
(TASEP) [7], asymmetric inclusion process (ASIP) [HI El HD] and many variations on 
these two models In all these models, particles jump between sites of a 

one- or higher-dimensional lattice and the dynamics is defined by specifying the hopping 
rates of the particles. The hopping rates are usually chosen so that there is a non-zero, 
macroscopic current of particles through the system driving it far from equilibrium, 
although the system often exhibits a non-equilibrium steady state independent of the 
initial condition. 

In this paper, we study an extension of the zero-range process to nearest-neighbour 
interactions, similar to that of Refs. [IS HZ]. In this model, particles interact when 
they are at the same site or at neighbouring sites. Although the model can be 
driven far from equilibrium, it is closely related to the equilibrium solid-on-solid (SOS) 
model [IS US [20] • A remarkable feature of this stochastic process is that the steady 
state factorises over pairs of neighbouring sites, also in dimensions higher than one, and 
thus we call it the pair-factorised steady state process (PFSS). This property facilitates 
analytical calculations in the one-dimensional version of the model, and in certain cases 
also in more than one dimension 1211 

Contrary to previous works which focused on generic properties of this model 
such as the existence of condensation [H] [22], the shape of the condensate HTJ , or 
generalisation to more complicated graphs 125 , we revisit here the original foundation 
of this model coming from solid-state physics, and choose a hopping rate which leads 
to the emergence of clusters of particles similar to the extended atomic islands known 
from non-equilibrium nanostructure formation [2S] and epitaxial thin film growth [22]- 
In these processes, a film of atoms is deposited on a substrate that serves as a template. 
One of three generic modes of epitaxial thin film growth [2S] - Stranski-Krastanov 
growth - has attracted considerable attention as it can be used, for example, to produce 
quantum dots [261 EZ]. In Stranski-Krastanov growth, deposited atoms form initially 
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Figure 1. Schematic stages of Stranski-Krastanov growth. Adatoms (blue spheres) 
are deposited on the substrate (purple spheres) until a desired density is reached, (a) 
Low density - an incomplete monolayer is formed, (b) As the density of adatoms 
increases, the adatoms form a complete monolayer and a partially filled second layer, 
(c) Upon further density increase, islands of variable height begin to form on the 
monolayer. Here the critical density of adatoms above which islands occur equals one 
adatom per one substrate site. The critical thickness depends on the mismatch between 
the substrate and adsorbate lattices, which is, however, not modelled explicitly in this 
paper. 


a flat, 2d layer. As the density of atoms on the snbstrate increases beyond a certain 
critical thickness, atomic islands start to nucleate as shown schematically in Fig. 

In this paper, we propose a simple, analytically tractable non-equilibrium toy 
model that mimics the 2d-to-3d transition observed in Stranski-Krastanov growth. In 
our model, we do not aim at reproducing all details of thin film growth (e.g., there 
is no mismatch between the substrate and adsorbate lattices) but we rather explore 
generic mechanisms that lead to island formation in non-equilibrium systems that mimic 
those of Stranski-Krastanov mode of growth. In particular, we show that by changing 
the strength of interactions between particles one can obtain different island shapes, 
similarly to what is seen in experiments. We also show that the shape is quite robust 
and does not change much when the system is pushed far from equilibrium either by 
imposing a macroscopic current of particles in one direction (as in electromigration on 
surfaces [28j), or by adding new particles to the system at a constant rate (the latter 
process imitating molecular beam epitaxy m)- 

2. The model 

The model that we study in this work comprises a two-dimensional, regular lattice with 
N = L X L sites and periodic boundary conditions in both directions. Let {m*} be the 
number of particles at sites i = 1,... ,N. The particles can be viewed as “adatoms” 
attached to the surface of the substrate where the number ruj corresponds to the height 
(in the third dimension) of a stack of atoms at the site i. We first consider the case 
when M particles have been deposited on the substrate and no further particles are 
being added, thus the number of particles is constant and equal to M. We shall later 
relax this assumption. 
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To model the dynamics of particles due to thermal excitations and external driving 
(e.g. electromigration), we assume the following rate at which a particle jumps out of 
site i: 

^ TT 9 {m - 

where {i,j) denotes all four nearest neighbours of site i and g{m,n) is a symmetric 
non-negative function to be specihed later. The jump rate depends on the number of 
particles at i and at its nearest neighbours, and by a suitable choice of g{m, n) we can 
replicate interactions between particles at neighbouring sites. The particle then hops 
to one of the neighbours with probabilities {r^} for k = 1 (right), 2 (left), 3 (top), and 
4 (bottom). The above choice of Uj is dictated by the requirement that the steady-state 
microstate probability assumes the following factorized form |22j . 


P(mi,... ,miv) 


^'[lg{mi,mj)d 

(ti) 


■ N 

mi — M 

_ ^=1 


( 2 ) 


where the probability factorises over pairs of neighbouring sites, the Kronecker delta (5[fc] 
(equal to 1 if A: = 0) ensures that the total number of particles is conserved, and Z is a 
normalisation constant giving §1 a valid probabilistic interpretation. The factorisation 
allows us to analyse the statics of the model using standard tools of statistical mechanics. 
Identifying P(mi,... ,mjv) with the Boltzmann distribution (1/2’) exp(—/5P) with the 
inverse temperature /5 = 1, we obtain the energy of the microstate 

P(mi,...,mjv) =-y^ln^(mi,mj). (3) 

(ti) 


Even though the system is in general out of equilibrium, the steady state is independent 
of the jump probabilities {r^}, and many steady-state quantities can be calculated as 
if the system was in equilibrium, with the microstate probability given by Eq. (§. 
The choice of the probabilities {r^} determines how far the system is from equilibrium; 
for example, for {rfc} = (1/3, 0,1/3,1/3} particles jump asymmetrically from left to 
right, which generates a macroscopic current of particles in this direction, whereas for 
{rk} = {1/4,1/4,1/4,1/4} the jumps are fully symmetric, the net current of particles 
is zero, and the system is at equilibrium. 

We also consider a (l-|-l)d counterpart of this model, in which particles jump to 
the right or left on a one-dimensional substrate, and the simplified form of the hopping 
rate Eq. Q is 




g{mi - 1, mj-i) g{mi - 1, rra^+i) 

g{mi,mi-i) g{mi,mi+i) 


The corresponding microstate probability then reads 


(4) 


P(mi,... ,miv) 


N 


^g{mi,mi+i)5 




■ N 

mi — M 

_ ^=1 


(5) 
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Similarly to the (2+l)d model, the probability ri = r 2 = 1/2 corresponds to the system 
in thermal equilibrium, whereas for ri = l,r 2 = 0 the particles can jump only to the 
right as in Ref. [Tb] . 

The model described above has been studied for a number of choices of g{m,n) 
in one dimension [ElIII!, and less extensively in two dimensions [21] and, given that 
g{m,n) fulfils certain criteria, the model is known to have a phase transition between 
a liquid and a condensed state as the density of particles crosses a threshold density. 
Depending on the choice of g{m,n), the condensate can be either localised at a single 
site as in the ZRP, or can be spatially extended over many sites nziEnj. 

In this paper, we aim to reproduce qualitatively the phenomenology of surface 
growth. We therefore assume the energy of a configuration (mi,..., miv) to be 


E 



( 6 ) 


which is equivalent to the following two-point symmetric weight function g[m,n): 


g{m, n) = exp 

where V (m) is the potential 

V{m) = U 


—J|m — n\ — -{V{m) + V{n)) 


X 9 / \ 3 

a \ I o 


m-|-l/ V?^ + l 


(7) 

( 8 ) 


The term proportional to J in the above expression represents the energy cost of “broken 
bonds” between neighbouring, vertical stacks of adatoms, whereas the term proportional 
to E accounts for interactions between adatoms and the substrate, see Fig.|^ We assume 
the latter to be described by the Lennard-Jones (LJ) potential, with a unity added to 
the denominator to make the expression finite for m = 0. Equation (|^ is relevant to 
interactions of molecules with a crystalline surface [30l |3T] , and the exponents 9 and 3 
arise from integrating the 12-6 LJ potential over the substrate surfac^ 

The model has three parameters that are related to adatom-substrate interactions: 
U,J and a. Large J suppresses, through |mj — mf, rapid variations in the height of 
neighbouring stacks {i,j) of adatoms, and flattens out the surface; large U makes the 
adatoms bind stronger to the substrate; a has the interpretation of the interaction range 
between the adatoms and the surface, measured in the units of the lattice constant. 
Figure shows how the potential V (m) behaves for different values of the parameter 
cr, for U = 1/2. Although our model only serves as heuristic means, the values of a 
we use throughout the paper (0 < cr < 3) fall in the range typically encountered in 
Stranski-Krastanov growth, see [Appendix A , 


:|: The model can be easily extended to other choices of the potential, including other exponents, and 
will exhibit qualitatively the same behaviour. 
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AE=-2J 



AE=V(3)+V(1)-2V(2) 



(V(3P\ 



Figure 2. When a particle jumps to a neighbouring site with rate the energy 
of the system changes. Left: a particle merging with a two-particle island reduces the 
number of height differences (“broken bonds”) — mi+i| by two; right: a particle 

hopping onto a neighbouring stack of adatoms changes the value of the on-site potential 
V{m). For cr = l,J>0, >0 both depicted energy changes are negative, AE < 0, 

hence these transitions are more likely than the moves in the opposite direction. 
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Figure 3. (a) Examples of the on-site potential V{m) of Eq. ([^ for {7 = 1/2. (b) 

Approximate potentials for cr = 0.5,1, and 2.5 constructed from one (dashed lines) or 
two (continuous line) delta functions; see Sec. 4.1-C. The plots have been vertically 


shifted for clarity; black horizontal lines mark the zero energy level in each case. 


3. Numerical results for the (2+l)d model 

In this section we discuss steady-state properties of the (2+l)d model with hxed number 
of particles. Since all quantities discussed here depend only on g{m,n) through the 
steady-state probability (j^ and not on transition probabilities between the states of 
the system, we took the liberty of using a Monte Carlo algorithm (for more details, 
to simulate the (2-|-l)d model on a computer. This approach, 
unsuitable for dynamic quantities such as the average current, is much better suited 
for simulations of large systems due to its significant speed gain over the dynamics 
described by Eq. 0 . 

Figure shows steady-state snapshots of the system for different surface densities 
p = M/L^, fixed U, J, and for two different a = 1,3. For cr = 1, increasing the density 
p produces first a fiat, irregular droplet of height m = 1 and size increasing with p; 


see 


Appendix B 
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Figure 4. (Color online) Snapshots of the steady-state of the system for different 
densities of particles p = Mll?\ colours represent different numbers of particles per 
site (see bar on the right). Each square represents a 64 x 64 lattice with periodic 
boundary conditions. The condensate (a green/yellow/red shape) forms for large 
enough densities. The parameters are J = !.!,[/ = 3, and (left) a = 1, (right) 
(7 = 3. 


then, above a certain critical density a hemi-spherical island - which we shall call 
the condensate - begins to form on the surface. The height of the condensate increases 
with p, while the height of the surface remains constant and equal to one. The situation 
looks similar for a = 3 except that the condensate forms on a layer of three particles 
thick. The snapshots suggest that the critical density for condensation is approximately 
equal to a, the range of the LJ potential. Indeed, simulations made for [/, J, a as in 
Fig. a and for a range of densities p = 1,..., 7 show that Pc = 1.0135 ± 0.0013 for a = 1 
and Pc = 3.068 ± 0.023 for cr = 3 (see Appendix C). The values of pc remain very close 
to [aj (i.e., the floor of a), which indicates that, for U large enough, the critical density 
is very close to the density of particles necessary to populate the first [crj layers above 
the substrate. We shall see in Sec. that this is also true for the analytically solvable 
(l+l)d model. 


3.1. Phase diagram 

To explore the effect of J, 17, and a, we have made simulations for fixed density p while 
varying J, U, a. Figure shows a pictorial representation of the U, J-phase diagram, for 
cr = 1. The density of particles is p = 3. The upper right corner of the diagram (large 
positive U, J) corresponds to the parameter region where particles condense into islands, 
and the bottom left (small J and negative U) to the region where only a fluctuating 
“wetting layer” can be observed. At the transition region between the condensate and 
the wetting-layer phase vertical or horizontal “stripes” of particles can be seen. They 
are caused by periodic boundary conditions: for small 17, due to a larger extension of the 
condensate its opposite sides merge together and form a stripe; a larger lattice would 
prevent this finite-size effect and hence the region where the stripes occur in fact belongs 
to the condensed phase. Stripes are also present for sufficiently negative 17, but their 
origin is different: the surface now repels adatoms, which tend to cluster together. 

The situation is qualitatively similar for other values of a, see Fig. For cr > 1, 
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Figure 5. (Color online) Phase diagram for a = I; each square represents a snapshot 
of the (2-|-l)d system with 64 x 64 sites, simulated for a given pair U, J. Colours 
represent different numbers of particles per site (see the colour bar). The condensates 
have been shifted so that each appears in the centre of the lattice. 



Figure 6. (Color online) Phase diagrams for a = 0.5,0.83, and 3 (left to right). 
Colors represent different numbers of particles per site (see the color bar). 
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U=3, J=0.5 U=3, J=1.1 
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Figure 7. (Color online) Top row: condensates on a square lattice rcy-plane) take more 
rectangular shapes for larger values of J (tr = 1, L = 200 in both cases). The colour 
coding has been chosen so as to enhance visibility of droplet boundaries. Circles (solid 
lines below) correspond to half of the maximal height. Bottom: central cross-sections 
of the condensate (a: 2 :-plane) for the same J, U as pictures in the upper row. The profile 
of the condensate becomes more rectangular as J increases, the height decreases, and 
the width in the ai-direction increases (see Appendix B). 


condensation occurs for ?7 >0; the surface is covered by a “wetting layer” composed 
of [aj layers of particles, with the condensate on the top-most layer. For cr < 1, the 
condensate (which again happens for U > 0) is surrounded by empty sites and there 
is no wetting layer. This is due to the fact that the potential V{m) has a minimum 
at m = 0, and is not sufficiently deep for m > 0 (cf. Fig. |^a)) for the particles to be 
attracted enough to the substrate. Similarly, for a > 1 and U < 0, the surface is empty 
apart from a localised, very high condensate peak (white spots in Fig. right panel; 
for high J there are no peaks because large surface tension does not let the simulation 
leave the flat initial condition). 

Figures |5]-[^ show that the shape of the condensate in the xy-plane depends on 
the parameters J, U : large U makes the condensate narrower and higher, while large 
J makes its surface (in the ; 2 -direction) flatter. Figure top row, shows that on a 
square lattice, small J leads to circularly-shaped condensates in the xy-plane, whereas 
for large J the condensate assumes a more square shape, reflecting the symmetry 
of the underlying lattice. The same figure, bottom row, demonstrates that as the 
condensate becomes more rectangular, its profile (section through the centre in the 
xz-plane) changes from an approximately parabolic to a more rectangular one. We 
shall see in Sec. E3 that the x. 2 -proflle of the (2-|-l)-dimensional condensate can be well 
approximated by the (l-|-l)d model. 

The shape of the condensate does not depend on whether the system is in 
equilibrium (by making the hops symmetric: = 1/4), or not. For example, even in the 
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extreme case when particles are allowed to jump to the right and not to the left (which 
produces a strong current in the a:-direction), the steady-state shape remains unaffected. 
This is caused by the lack of any explicit dependence of the steady-state distribution 
Q on the hopping probabilities {r^}. The dynamics of condensation, however, will be 
different for different {rfc}. 

Before we discuss the dynamics of condensation, let us briefly comment on the 
relation between what we see in our model, and experiments on Stranski-Krastanov 
growth, which partly motivated this study. As explained in the introduction, Stranski- 
Krastanov growth can be used to produce quantum dots. A visual comparison between 
our results, and experimentally obtained Ala;Gai_a;As quantum dots on GaAs [22! 
(Figs. 9-11 therein) and GaN on AIN |27| (Fig. 2 therein) shows that these quantum dots 
are not dissimilar to our condensates, and their shape depends on growth conditions and 
chemical composition which corresponds to different values of J, U parameters in our 
model. Our model can therefore qualitatively reproduce certain aspects of the growth 
of quantum dots. 

3.2. Dynamics 

We now discuss the dynamics of condensation in this model. The equilibrium Monte 
Carlo algorithm employed in the previous section cannot be used here and we have 
to simulate the process using a kinetic Monte Carlo algorithm with the hopping rate 
(§. Figure left, shows the time evolution of the equilibrium model (r^ = 1/4), 
starting from randomly distributed particles at t = 0. The process has two phases. 
First, particles rapidly aggregate into clusters; the second, slower stage involves clusters 
exchanging particles through the background. Eventually, only one cluster - the 
condensate - remains in the system. 

The time it takes for a single condensate to build up can be found using similar 
arguments to those for the ZRP [3|. Numerical simulations suggest that the time 
T to condensation is dominated by the process of merging the last few remaining 
clusters. Each cluster has on average 0{M) = 0{{p — Pc)L‘^) particles and the inter¬ 
cluster distance is 0{L). A cluster of size m loses particles through its boundary. 
The rate Uemit with which each of the I = l{m) sites at the edge of the cluster 
emits particles only weakly depends on the size if m ^ 1. For example, the rate 
at which particles are emitted from a site of height h in the condensate’s wall is 
Uemit = {g{h-l,h)/g{h,h)f {g{h-l,pc)/gih,pc)) = exp(-2J) for any /i > 1, and 
hence we can take it to be constant for all clusters. The total emission rate of the 
cluster is /uemit- Put differently, each such particle is emitted every Temit = l/(^Uemit) 
time units. Once they leave the cluster, the particles undergo a random walk with 
diffusion constant D ^ u{pc + l\pc, Pc, Pc, Pc) = [diPc Pc)/9{Pc + l;Pc)]^- Most of these 
particles are quickly reabsorbed due to recurrence of 2d random walk |32| but particles 
that have departed a distance 0(L) can be intercepted by other clusters. The time 
the particle needs to travel to reach another cluster is approximately Ttravei = 0{L‘^/D). 
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Figure 8. Example of the time evolution of a 32 x 32 system, for J = 1.1, U = 15, cr = 
l,p = 2. Each xy plane (with periodic boundary conditions) corresponds to a single 
time frame. Cubes represent occupied sites with two or more particles (the background 
has 1 particle per site). Left; equilibrium model (r* = 1/4); right; non-equilibrium 
model (ri = rs = r 4 = 1/3, r 2 = 0). 


Since Ttravei increases with L whereas Temit does not, diffusion is the limiting step and the 
total time it takes to move a particle from one cluster to another one is approximately 
equal to Ttravei for large enough L. 

Large clusters e mi t more particles, but they also re-absorb them with higher 
probability due to their larger circumference. A particle that has diffused away is more 
likely to be absorbed by a large than a small cluster. This causes a net current of particles 
flowing from smaller to larger clusters. The time it takes to transfer 0{Lf) particles 
from a small cluster to the condensate is therefore T = O(T^Ttravei) = 0{L‘^/D). The 
scaling T oc D~^L^ is verified in Fig. left, where we plot the average time it takes to 
have only one cluster (the condensate) in the system. 

Interestingly, the non-equilibrium case in which there is a current of particles in the 
y direction, leads to the same prediction (Fig. right). Although individual particles 
drift preferably in the y direction, the clusters remain quasi-static, see Fig. right. 
Moreover, the time for a particle to move from one cluster to another is again O(L^); 
its motion is still diffusive in the x direction and, unless the clusters are accidentally 
aligned so that particles can move between them in straight lines, diffusion dominates 
over the ballistic motion for which the time scale would be 0{L). 

4. Analytically soluble (l+l)d model 

Although the steady-state probability of the (2-|-l)d model has a simple, factorised form 
([^, its exact solution remains elusive. However, we can learn about many properties 
of this model using approximate methods. In this section we shall analyse a (l+l)d 
counterpart of the model, described by Eqs. 0 ,@. We shall show that this simpler 
model predicts not only the critical density Pc but also qualitatively reproduces the 
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Figure 9. The time to condensation for different sizes L = 24,32,48,64,96. Left: 
equilibrium model {vk = 1/2), right: non-equilibrium model (ri = rs = r 4 = 1/3, r 2 = 
0). In all cases J = 1.1, = 15, a = 1. The solid line is T = D~^L^ with 

D = (^(1, l)/^(2,1))^ = 1076.88 for the left panel and D = 4 x 1076.88 for the right 
panel. Note that the unknown proportionality coefficient in the anticipated asymptotic 
formula T oc D~^L^ turns out to be close to one in the non-equilibrium case and to 
four in the equilibrium case. 


transition lines in the {U, J)-phase diagrams from Figs. 

The model can be analysed along the same lines as in Ref. [22]. Let us define the 
canonical partition function 


L 

Z{L,M) = ER g{mi,mi+i)S 

{rrii} i=l 


■ L 

mi — M 


(9) 


and its grand-canonical counterpart 


L 

Zl{z) = '^Z{L,M)z^ = zZ^i^iY[9{mi,mi+i), ( 10 ) 

M {TUi} i=l 


where .2 is the fugacity, determined from its relationship to the average density, 


p{z) = 



2 ; dlnZL^z) 
L ^ 


( 11 ) 


Thanks to the one-dimensionality of the problem, Zl{z) can be expressed using the 
standard transfer-matrix approach: 


Z, 


.w= E 


T T 

mim2 m2m3 


Tm^mi = TyT{zY 




( 12 ) 


where T^n = denotes the transfer matrix. We expect the partition 

function Zl{z) to have a finite radius of convergence 2c. If p{z) 00 as z increases 
from 0 to 2c, the grand-canonical ensemble is valid for any density of particles, and the 
probability p{m) of finding m particles at a randomly chosen site reads 


p{m) 



,2 ’ 

Z-^m=0 


(13) 
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where is the principal eigenvector (associated with the largest eigenvalue) of T^niz)- 
If, however, p{z) has a finite limiting value p{z) — p^ as 2 approaches z^, then Eq. ( [TT| ) 
cannot be satisfied for z > Zc (or equivalently for p > p^) and the grand-canonical 
ensemble cannot be constructed. This corresponds to a transition from the liquid to the 
condensed state for p > p^. 

Since for our choice ([^ of g{m,n) the critical Zc = 1, to determine the critical 
density at which this transition happens, we must find the eigenvector fra of the matrix 
Tm,n = n) to the maximal eigenvalue Amax- 

g(m, Tl)(j)ri ^maxfrm 


and, rewriting the partition function in the large-L limit as Zl{z) = from Eq. (0 

we obtain 


Pc = ( 15 ) 

Em=0 ^ra 

Condensation can occur only if fra decays with m faster than ~ otherwise the 

critical density pc becomes infinite in the thermodynamic limit. It turns out that 
although the eigenproblem (14) is very easy to solve numerically, it is still too difficult to 
solve analytically for our particular choice of g{m, n) from Eq. 0 . To make progress, 
we observe that since the occupation numbers are discrete, the potential values are 
discretised as well. Moreover, the value of V (m) varies significantly with m only for the 
first few integer m. This suggests that the potential can be approximated by a sum of 
a few Kronecker delta functions with appropriate amplitudes. 


4 . 1 . Solution for a <C 1 

Let us first consider the case a <C 1. As illustrated in Eig. [^a), the on-site potential 
V{m) has one dominant minimum. This allows us to approximate the potential as 
— US[m], where 5[m\ is the Kronecker delta, and U = U{a^ — a^). The value of IJ is 
chosen to reproduce the value obtained from the exact formula (|^ for m = 0. Eigure|^b) 
shows this approximate potential for a = 0.5. The term 5 [m] lowers the energy and 
hence it increases the probability of a state in which the occupation is m = 0; non-zero 
occupation, on the other hand, is energetically unfavourable. Physically, this could mean 
that the particles cannot wet the substrate that is strongly “hydrophobic”. Therefore, 
the on-site potential favours empty sites, which leads to mass condensation seen as an 
“island” of particles discussed before. 

This model, which we shall call “model A” here, can be solved using the same 
approach as in Ref. [22] • Assuming the weight matrix 

g(rn, n) = exp[—J|m — n\-\-U{5 [m] -I- 8 [n])/2], (16) 


the eigenvector fra which solves Eq. 0 must take the form fm oc exp(A5 [m] -|- Bm), 
with some constants A,B. Inserting it into Eq. (14) we obtain the constants A = 
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U/2,B = —J — ln[l — exp(—t/)] = —J + Jq, where Jo = — ln[l — exp(—t/)]. Since 
(frn exp(5m), the entries of increase with increasing m for J < Jq and, from 


Eq. (15), the critical density pc is infinite. This means that condensation cannot occnr 
even for a very high density p of particles if J < Jq. For J > Jq, however, <prn falls off 
exponentially, and the critical density for model A reads 


Pc = 




,Jo _ I 


TZ=o C _ i]' 


(17) 


For example, for a = 0.5 and U = J = 2 (where J > Jq ^ 1.5), the critical density 
Pc ~ 0.53. The transition line, which separates the region in the {U, J)-plane where 
condensation occurs from the region where it does not, is given by 


J = — ln(l 






(18) 


In Fig. 10 we show a plot of Eq. (18) compared with an exact, numerical solution to the 


eigenproblem (14). The agreement is good even for a relatively large a = 0.5. It is also 


worth noticing that the transition line predicted by this model is qualitatively similar 
to that of the (2+l)d model from Fig. 


4-2. Solution for a ^ 1 

In the case a ^ 1, the minimum of the potential is located at m = \_a\. Since for all 
occupation numbers m < [aj the potential is very large (see Fig.|^, we can approximate 
it by 

—U5 [m — [aj], for m > [crj , 
oo, for m < \_a\ , 


V{m) 


(19) 


where 


U = U 


a 


L^J +1 


a 


[<tJ + 1 


( 20 ) 


Let us call this “model B”. The potential barrier at mo — [crJ — 1 means that if the 
density of particles p > [aj, no site will have less than [crj particles. Consequently, 
all elements m = 0,... , [crj — 1 of the eigenvector will be zero. We can derive the 
critical density for condensation in the same way as in the previous section. Assuming 
that (fra oc exp(A(5 [m — [crj] + Sm) for m > [crj , we obtain that A = U /2, B = —J + Jq, 
where Jq = — ln[l — exp(—C)] has the same form as previously, albeit with a different 
U given by Eq. (20). The critical density given by Eq. (17) is shifted by [crj which 
accounts for the vanishing elements of the eigenvector (the shift by x occurs for any 
g(m,n) whose eigenvector’s first x elements vanish): 


Pc = + Pc = + 


;Jo _ 1 


— g—2(J—Jo)j |^g2(J—Jo) — ]_] ' 


( 21 ) 
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Figure 10. (Colour online) Phase transition lines for the (l+l)d model with 
LJ potential The dashed lines are analytic solutions: The long-dashed curve 

corresponds to model A [Eq. (18)], the short-dashed and the dot-dashed to the double¬ 
delta model C [Eq. ^7\]. The points are numerical solutions obtained by diagonalising 


g{mm) numerically (see Appendix C). 


The critical line J{U) separating the condensed and liquid phases is 

3 / \ 9^ 


J = — In < 1 — exp 


-U 


a 


L(tJ +1 


a 


kj+ 1 , 


( 22 ) 


This expression is a good approximation for the critical line for large a, but is much 
worse for a ~ 2 — 3 that we use here in simulations. In the next section, |4.31 we show 
that if the potential V{m) is approximated by a sum of two delta functions (“model C”), 


the agreement between the approximate solution J{U) [see Eq. (27)] and the simulation 
data becomes much better for smaller a, as seen in Fig. [iq 


4 . 3 . Double-delta approximation of the potential 


In previous sections we used a single-Kronecker delta approximation of the LJ potential. 
This approximation is fairly sufficient for the two limits from Secs. igB for (T -C 1 and 
cr ^ 1, but it does not work well for a ~ 1 which is the range we are interested in in 
this work. The approximation can be improved by modelling the potential with two 
Kronecker delta functions with suitable amplitudes: 


Vim) 


—UiS [m — [aj] — U 2 S [m — [aj — 1], for m > [aj 
00 , for m < [_a\ 


(23) 


where 
Ui = U 


[ctJ + 1 


[(tJ + 1 


and Uo = U 


[crj -S 2 


a 


[(tJ + 2 


( 24 ) 
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The potential (23) is shown in Fig. |^b) for a = 1 and 2.5 (blue continuous lines). 
We shall refer to this model as “model C”. The principal eigenvector of g(m,n) 
for this model is given by oc exp(74i5 [n — [crj] + ^42^ [n — [crj — 1] + Bn) (and 

(po = ■ ■ ■ = = 0), where ^42 = b^2/2. The parameters Ai and B can be determined 

analytically from the third order polynomial equation in 


Q ^ ^At+Utl2 gB+J+(72 _ 


e2(s+J) 

- 1 ’ 


(25) 


on substitution of Ai 

Ai = ^ + 2B-In [e^+^ + (l - 2e^ cosh J + e^^)] . (26) 

The transition line, which corresponds to the condition B = 0, can be determined from 
a cubic equation in e'^, 


('i _ ^J) ^ gJj (^2-2 cosh J + . (27) 

This equation can be solved exactly, but the formula for J = J{U) is complicated and 
not very illuminating, hence we omit it here and only plot the solution in Fig. [T^ 


f.f. Shape of the condensate 

Above Pc, a spatially extended condensate forms in the system, see Fig. The figure 
shows that the shape of the condensate, obtained from MC simulations by shifting the 
condensate to i = L/2 and averaging over many samples, is approximately parabolic. 
The shape can be analytically derived using the result of Ref. [22] for a (l+l)d model 
with a weight function g{m,n) = K{\m — n\)^p{m)p{n), where K{m) and p{m) are 
arbitrary functions decaying sufficiently fast with m ^ oo. Let us assume that the 
condensate has mass M' = M — pcL, and define rescaled variables h = {mi)/'/M' 
and t = — 1, where (m*) denotes the mean occupation of lattice site i and w 

is a constant which we shall determine later. In the large-L limit, the shape of the 
condensate in these new variables is given by [22] 


h{t) 


w Kiv) 

TT- ffi “=-) 

2v K{vf) 


(28) 


where K{x) = , and w and v are auxiliary parameters that can be 

obtained from A^ax = K{v) and 


W = V 



xk'{x) 

k{x) 




(29) 


The one-point weight function p{m) enters these formulas only through the largest 
eigenvalue A^ax of the matrix g(m,n). Applying these results to our model with LJ 
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Figure 11. The shape of the (l+l)d condensate plotted in the normalised variables 
(t,h) for cr = 1. The curves represent Eq. (30), where for black continuous lines Amax 
was obtained by numerical diagonalisation of a 100 x 100 ^(m, n) matrix, and for blue 
dashed line by solving model C ( the two curves are identical in the right panel); the 
circles come from a MC simulation of LJ system of size L = 2000 with M = 60000 
particles, averaged over 10^ MC sweeps. The layer of thickness pc ~ 1.6 and 1.0 
has been subtracted in the left and right panel, respectively. The actual height and 
width of the condensate are h{0)\/M' « 65 and w\/M' 1302 for the left panel, and 

h(0) VM' 48 and « 1379 for the right one, respectively. 


potential, we have K{m) = e and hence the function K[x) 
shape h{t) reads 

, , , w ^ (cosh J — cosh vt \ 

hit) = — In ---^— , 

2v \ cosh J — cosh v J 

where v must be determined from the equation 


sinh J 

cosh J—cosh X ' 


The 


(30) 


A 


max 


sinh J 

cosh J — cosh V 


(31) 


As already mentioned, the only dependence on the potential V{m) = — Inp(m) is 
through the eigenvalue Amax of g{m,n), which can be found numerically for the LJ 
potential, and analytically for the approximate models A-C, for which 


V = J - Jo, 


(32) 


and Jo is obtained from Eqs. (18), (22), and (27), for the respective models. Equation 


(30) is a good approximation to the exact shape of the (l+l)d condensate already for 
relatively small M, see Eig. [m 

It is also interesting to note that the two-dimensional sections of the (2-|-l)d 


condensate resemble very closely the (l+l)d envelope. In Fig. 12 we show that the 
a:. 2 -section through the centre of the condensate as well as further non-central sections 


are very well approximated by Eq. (30), with v, J htted to the numerical data. However, 


we do not know whether this similarity is not a mere coincidence, nor could we find an 
analytical formula which would predict the “effective” constants w, v from the “bare” 
parameters U, J of the (2-|-l)d model. 
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Figure 12. Cross-sections of the (2-|-l)d condensate with LJ potential {a = 1) 
in xz plane (circles). The biggest envelope corresponds to the most-central section 
[y = 0), with other sections taken at y = 10,20,30,... lattice sites. In sum, 8 
sections are plotted in the left and 9 in the right panel. The continuous lines are the 
(l-|-l)d theoretical shapes (30) fitted with only two parameters v, J, and w = w{v, J) 
obtained from Eq. (29). The parameters v,J were fitted separately for each section; 
the “effective” parameters J and U = — In [(1 -|- 1)/— 1)] (an exact solution for 
the (l-|-l)d delta model) obtained from the fit decrease monotonously with increasing 
distance y from the central section. The details of the simulation can be found 
in [Appendix B| 


5. Deposition of new particles with constant rate 

One of the features of our non-equilibrium model is that its steady-state probability 
assumes a relatively simple, factorised form ([^ and, as we have seen, this allows us to 
calculate some quantities analytically. In this section, we explore the consequences of 
breaking this factorisation by releasing the constraint of mass conservation. 

In the new model, particles are added to the system at a constant rate a, as in 
molecular beam epitaxy. This model does not have a steady state in the sense of the 
constant-mass model from previous sections, because the number of particles per site 
increases over time. However, we shall see that the model has a quasi-steady state when 
the number of deposited particles is not too large, and that this state is very similar to 
what we discussed before. 

Figure shows snapshots of the system at different times, for two different (low 
and high) mass deposition rates. For low deposition rates it can be seen that a single 
condensate is formed. This is not unexpected: particles jump between lattice sites much 
faster than it takes to add a new particle, and the system relaxes to a quasi-steady state 
similar to that of the constant-mass model. However, when the deposition rate a is 
high enough, new condensates are formed faster than they can coalesce. In this regime 
multiple condensates arise. 

We can estimate the magnitude of the deposition rate ttgep that separates the two 
regimes as follows. We consider only what happens after first [aj layers have been 
filled because this is when condensation begins. A newly added particle stays on the 
surface of the top-most layer and performs a random walk with diffusion constant D 
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(see Sec. 3.2) until it collides with another particle and becomes the seed of a new 
cluster. Let us denote (with a slight abuse of notation) the quasi-steady state density of 
such isolated particles by p[§} This excludes particles from the complete layers as well 
as particles in the clnsters. If we neglect spatial correlations, the probability that onr 
particle collides with another one during the next step is p for p 1. The probability 
that the particle has not yet collided after n steps is then (1 — p)^, and the mean 
number of steps to collision (n) = '^j^np{l — pY = Y ~ p)/p — Vp- The time to 
collision is then ~ l/(Dp). Dnring this time the particle departs from its starting point 
by (r) ~ a/( n) ~ p This distance gives the characteristic length scale for spatial 
separation of clusters of particles. If it is of the order of the spatial extension of the 
simulation box L, only one cluster — the condensate — will form in the system. By 
equating (r) and L we obtain the density p ~ 1/L^ at which this happens. To relate 
this density to the deposition rate a we note that deposition mnst be balanced by 
the rate with which particles form clusters; since the clusters are relatively narrow, 
their contribution to the average density of particles can be neglected. This gives 
us a = pL^Dp where pLf is the number of “free” particles in the system. Inserting 
p ~ 1/L^ we obtain ttgep D/L'^. Hence, if agep D/Lf, multiple condensates are 
present in the system, otherwise there is only one condensate. Figure IT shows the 
inverse participation ratio (IPR) of the occupation numbers {ruj}, which approximately 
corresponds to the nnmber of condensates, as a fnnction of the density p of the already 
deposited mass (proportional to the physical time), for different deposition rates a. 
The figure indicates that the theoretically predicted Osep correctly estimates the critical 
deposition rate if ttgep ~ 5D/L^, i.e. the proportionality factor is of the order of 5. 

Regardless of whether the deposition rate is high or low, the shape of the 
condensate(s) can still be well approximated by the equilibrium (l-|-l)d analytical 
solution. This is illustrated in Fig. [T^ where we compare the shape obtained in 
simulations of the (2-|-l)d LJ model for mass deposition rate a — 0.22 for J — 0.5 
and a = 0.62 for J = 1.1 to the exact solution (30) for the (l+l)d model, with w,v 
fitted to the cross-sections of the (2-|-l)d condensate (see Appendix B for more details). 


6. Conclusions 

In this work, motivated by thin-film growth processes and, in particular, by the 
Stranski-Krastanov growth mode, we propose a simple, non-equilibrium physics model 
in which spatially extended condensates (“islands”) form when the density of particles 
exceeds a critical value. Our model assumes short-range, valence-bond type interactions 
between particles, and Lennard-Jones interactions between particles and the substrate 
on which the growth occurs. Depending on the range a of the Lennard-Jones potential, 
condensation occurs either directly on the substrate (for a < 1) or on a previously 
formed layer of several particles thick (for a > 1). 

§ This new p should not be confused with p = M/LP defined previously. 
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Figure 13. (Colour online) Simulations of a non-equilibrium system with J = 1.1, = 

3, cr = 1 on a 128 x 128 lattice for a constant rate of mass deposition: (top row) a = 9.85 
incoming particles per unit time, (bottom row) a = 0.31 particles per unit time. The 
time ranges (top) from 1.5 x 10^, 1.8 x 10^, 3.3 x 10^ to 8.3 x 10^ and (bottom) from 
47.9 X 10^,58.5 x 10^, 106.4 x 10^ to 266 x 10^ time units. Multiple condensates form 
when the deposition rate is high enough (top), whereas for low a (bottom) only one 
condensate is created. 



Figure 14. (Colour online) The inverse participation ratio, IPR = (rfiffi j m? 
(where the sums run over all existing island masses), of condensates formed in 
J = 1.1, = 3,cr = 1 systems with different deposition rates (the fastest is the red 

topmost curve, the slowest is the blue bottommost one) as a function of the density p 
of the already deposited mass. The mass influxes are given in particles per unit time. 
The error bars are standard deviations of 20 simulation runs; the ^-axis is logarithmic. 
The theoretical estimate asep (x D/L^ for the rate separating regimes with one and 
many condensates yields 0.13 and 0.033 for L = 32 and 64, respectively. Assuming the 
proportionality factor 5, the estimated critical densities are 0.65 and 0.165, respectively, 
and they can be seen to separate well the curves for which the IPR remains very close 
to 1 for all densities (times) and for which it is larger than 1. 
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Figure 15. Sections of (2+l)d condensates as in Fig. 12, but with a constant 


influx of mass (0.22 and 0.62 particles per unit time for J = 0.5 and 1.1, where 
Usep ~ = 7.6 X 10-^ and 8.3 X 10 respectively). The circles are sections 

of condensate snapshots averaged over time and over the multiple condensates, with 
masses between about 12 - 34000 particles, that formed simultaneously in the system. 
The red dashed lines correspond to the fitted Eq. (30), whereas black continuous lines 
to the analytical shape corrected for the width fluctuations |22] responsible for the tails 
on the brim of the condensate. For details, see [Appendix B| 


Although there have been numerous approaches to simulating thin-film growth, 
(see, e.g., the review |H3], or a recent kinetic Monte Carlo study |H1!), the most 
interesting featnre of onr oversimplified model is that it enables us to calculate many 
quantities analytically. This is possible due to a pair-factorised steady state (PFSS) 
probability of microstates in our model. In a (l-|-l)d version of the model, we have 
been able to derive the phase diagram of the model, to calculate the critical density 
for condensation, and to find the shape of the condensate which turned out to depend 
on the strength of adatom-adatom and adatom-substrate interactions. In the (2-|-l)d 
model, which corresponds to the physically relevant growth of 2d layers of adatoms, 
we have shown that the shape of the condensate is well approximated by the (l-|-l)d 
solution. 

We have also studied an open system in which new particles are added at a constant 
rate. We have shown that condensation occnrs above a certain density of particles, and 
although it is a transient phenomenon, the properties of the condensate are similar to 
those of the model with mass conservation. 

In this work, we have focused on the steady-state, or quasi-steady state properties 
of the condensate and its late-time dynamics. It wonld be interesting to broaden onr 
research into the kinetics of initial steps of condensate formation. Further research could 
also involve manipulating the geometry of the underlying lattice, e.g., introduction of 
lattice defects which could imitate heteroepitaxial growth more closely [35l |36] . 
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Appendix A. Typical values of a for thin-layer growth 


Assuming that the substrate atoms are uniformly distributed over the lower half-space 
{z < 0) of the system, the form of the integrated LJ potential is [3T] 


27r na'^ M 

3 d 15 \ r/ V^/ ’ 


(A.l) 


where n is the number density of substrate atoms on the surface, e has the dimension 
of energy per mol, a' is the range of the LJ potential, and d is the layer spacing of the 
substrate. The parameters U and a from our formula (j^ can be expressed through 
n, d, o' as: 



Carbon or silicon crystals are usually modelled with, very roughly, 3A< o' < 
4A |37l |38l |39], whereas lattice constants of C, Si or GaAs are respectively 3.56A, 
5 . 43 A, and 5 . 65 A [301, which yields U ~ 0.05 eV and a ~ 2.5A. Taking into account 
that in our work all distances are measured in terms of the lattice spacing d, and that 
d ~ lA in most metals, o should be about 0.5 — 3; these are the values we use in this 
work. The value of J, on the other hand, can be approximated by the Ehrlich-Schwoebel 
barrier energy, which is typically of the order 0.1 — 0.5 eV. Together with k^T set to 1 
in our simulations and a liquid nitrogen cooled molecular-beam epitaxy temperature of 
77 K, we get the very rough estimates of U ps 10, J ~ 50. In our model, however, we 
use J ~ 1 because for significantly greater values the acceptance rate in MC simulations 
would become many orders of magnitude smaller, and consequently the simulation times 
would become unfeasible. 
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To determine the phase diagrams in Figs. [5]-[^ we used equilibrium Monte Carlo 
simulations with Metropolis acceptance probability mi. A single move consisted of 
picking up a random site and, if it was occupied, moving a particle to another randomly 
chosen site anywhere in the system. In comparison to the stochastic simulation of the 
original dynamics of the model, this significantly reduced the computation time while 
preserving the stationary state j2T]. For each pair {U, J), the 64 x 64 system (with p = 3 
for a = 0.5,0.8,1 and p = 6 for a = 3) was simulated for 4 x 10^ sweep^and, prior to 
that, it was thermalised for 2 x 10^ sweeps. The strongly rectangular shape of the island 
for high U and J values is due to the geometry of the square lattice and is independent 
of the initial conditions. 

In Figs. and 1^, Monte Carlo simulations were performed on a lattice of size 
A = 200 X 200 with M = 25A = 10® particles and the Lennard-Jones on-site potential. 
Both cuboid and cylindrical initial condition were used and as we did not find any 
differences between them, we concluded that thermalisation was long enough to erase any 
trace of the initial configuration. The simulations took 8 x 10^ time steps (around four 
weeks of computer time), half of which was thermalisation, for the cuboid (150 x 150x44) 
initial condition, and 4 x 10^ time steps, 10% of which was thermalisation, for the 
cylindrical (diameter 140, height 60) initial condition. The final plots shown were 
obtained from the latter simulation. 

The simulations of the dynamics of the (2-|-l)d model, and the model with mass 
deposition were performed using a simplified, kinetic Monte Carlo algorithm. Each time 
step a random site was picked and, if it was non-empty, one of the nearest neighbours 
was chosen with probabilities {ri,r 2 ,ra,r 4 } for right, left, top, and bottom jumps, 
respectively. The particle was then moved between these two sites with probability 
u/uma.x where u is the rate from Eq. 0 and Umax was chosen to be larger than the largest 
possible hop rate for a given set of parameters. This procedure was repeated times. 
Finally, the physical time was incremented by dt = 1/umax- In the model with mass 
deposition, a new particle was added every dt/a steps. This algorithm, although very 
fast, differs slightly from genuine kinetic Monte Carlo algorithms such as the Gillespie 
algorithm H- However, we checked that both algorithms produce indistinguishable 
results when averaged over a sufficiently long time. Simulations for Figs. and [T^ were 
performed on a 128 x 128 square lattice, with one particle at a randomly chosen site as 
the initial configuration. For simulations in Fig. with 32 x 32 and 64 x 64 square 
lattice systems, we counted as condensates all clusters both occupying an area greater 
than 1 site and having a height greater than 1 particle. The time to condensation in 
Fig. was determined as the average time (20-100 simulations per data point) at which 
the number of clusters larger than p = M/N dropped to one for the first time. 

The histograms of the (2+l)d condensates in Fig. 15 were obtained from a single 
simulation run with mass influx a = 0.62 for J = 1.1, a = 0.22 for J = 0.5. 


II A “sweep” comprises L attempted moves, whereas in (2+l)d it corresponds to attempts. 
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The condensate heights were rescaled according to their masses, the discrete lattice 
occupations were (linearly) interpolated, and only then were the interpolations averaged 
prodncing the histograms. 

The simnlations of the (l+l)d systems with cr = 1 for Fig. [IT| were performed on 
L = 2000 nodes with M = 60000 particles. The simulations took 8 x 10^ sweeps, with 
7 X 10^ sweeps devoted to thermalisation, and 10^ for recording the histogram. The 
theoretical Pc ~ 1-581 and the actual subtracted background was p ~ 1.599 ± 0.004 
thick for U = 1.5, J = 2(pc = p = l± 0.001 for U = 4, J = 8). The theoretical height 
of the condensate was h{0)'/M' ~ 64.84, and the actual height measured in simulations 
was 65.03 {h{0)'/M' ~ 48.02, simulations: 47.38). 


Appendix C. Critical density 


The transition lines shown for the Lennard-Jones potential in Fig. 10 were obtained 


numerically by diagonalising the matrix g{m,n), as in Eq. (14). For faster performance, 
only a 21-element wide band was retained in the matrix (10 elements below and above 
the diagonal; the furthermost elements are of the order of exp(—lOJ)), but to avoid 
numerical errors we used a direct banded matrix solver instead of the iterative (e.g., 
Lanczos) method. The parameter U was sampled at 0.025 intervals and the parameter 
J was determined by the bisection method (the last step of size AJ = 0.0195). The 
points where the critical density pc from Eq. ii increased slower than a logarithm of 
the matrix size were considered to belong to the condensed phase. The behaviour was 
classified as either slower or faster than logarithmic by: first, measuring Pc(T) for the 
weight matrix sizes 250,500,1000,2000; next, fitting a line in InL for the first three 
points, and another one for the last three points; finally, comparing the two slopes and 
if the second one was lower, classifying a given U, J pair as belonging to the condensed 
phase. 

In order to determine the critical density of particles above which condensation 
occurs, we simulated the model with fixed J, U,a while varying the density p = M/L. 
Each simulation was thermalised prior to measuring the mass M' of the condensate. 

We then used linear regression M' = M — Lp^ to determine pc from the sizes of 
the condensate for different M’s, taking into account only sufficiently large M'’s, see 


Eig.[^ 

We also performed simulations close to the expected Pc, as shown in the insets of 
Eig.|^ The results indicate that there is a non-linear drop in the condensate mass 
near pc and hence our method may have produced small but systematic errors when 
estimating the critical density via linear regression. 
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